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Abstract — The Brazilian productive sector of pulp and 
paper represents a relevant contribution for the 
development of Brazil. To increase the competitiveness of 
Brazilian companies to an international level, products 
must have high standards of quality and high added value. 
Thus, the mathematical modeling and simulation of 
industrial processes ensures the stability of production. 
This study presents the fit of mathematical models for the 
Oxygen Delignification process of eucalyptus pulp of the 
industry Klabin Monte Alegre. The mathematical model 
estimates the kappa number after the reactor, based on 
two kinetic models given by the literature, one of these 
models considers oxygen excess in the reaction medium. 
The models showed a mean relative error of 10 %. The 
adjustment of the kinetic parameters equations was done 
in Matlab software, using classical methods of 
optimization, such as BFGS, DFP, Steepest Descent, 
Gauss Newton, Simplex and Levenberg Marquardt. The 
models were incorporated in the commercial simulator 
CADSIM Plus to provide an optimization tool to the pulp 
industries. The simulator predicts the kappa number after 
the Oxygen Delignification reactor. The results of the 
phenomenological models indicate that possibly there is 
excess of oxygen in the reaction media. Only the model 
that considered the presence of the oxygen in the kinetic 
equation was able to be implemented in the simulator 
CADSIM Plus, generating consistent results, with an 
absolute error of ±2 kappa number. 


Application: The kinetic model applied to the CADSIM 
Plus software in this study may be used to optimize the 
Oxygen Delignification process either by reducing 
chemical consumptions or by testing different process 
conditions without changing production. 

Keywords — Cellulosic Pulp, Kinetic Expressions, 
CADSIM Plus, 

I. INTRODUCTION 

Klabin Papeis Monte Alegre increased the paper 
production capacity of the mill to 1,1 million tons/y. This 
project gave the company an excellent position in the 
global market with highly competitive production costs. 
The incorporation of mathematical models and process 
simulators helps the mill to optimize the pulp production 
assuring quality specifications that fit the clients demand 
with a safety and environmentally friend process. 

After the Kraft cooking in the continuous digester 
# 1, the pulp goes to the oxygen delignification process 
that occurs in a medium consistency media with sodium 
hydroxide as alkali source. The pulp is washed and 
discharged into a tank where alkali and oxygen with 
medium pressure steam are applied. There are two parallel 
lines with two post oxygen washers in each line. At one 
line, the pulp is washed in a pressurized diffuser washer 
followed by a post oxygen filter washer. In the other line 
the washing is done with two presses. The washed pulp 
goes to the ECF bleaching. Figure 1 illustrates the mill 
process. 
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Fig.l: Oxygen delignification plant of Klabin Monte Alegre. 


The oxygen delignification is the process which 
the pulp, suspended into an alkaline solution, is 
pressurized with oxygen which forms a stable dispersion in 
the pulp and is consumed in the lignin reactions. This stage 
removes 50% of the reminiscent lignin in the brown pulp, 
saving chemicals in the Bleach Plant and minimizing 
organic material in the effluents (COD and BOD). The 
process is controlled by a heterogeneous reaction that 
involves the solid fibers, liquid phase around and between 
the fibers and oxygen dispersed gas phase [1]. In general, 
delignification and the carbohydrates degradation reactions 
are mainly affected by temperature, sodium hydroxide 
concentration and dissolved oxygen concentration. 

Rubini [2] developed mathematical models based 
on the phenomenological and neural network models to 
predict the kappa number after the oxygen delignification 
reactor. The results achieved were very satisfactory 
regarding industrial simulation having an average error of 
8,5%. By those results and the similarity of the data used 
in this study, the same kinetic models used at Rubini’s [2] 
work were used given by Agarwal et al [3] and Violette 
[1]. The model proposed by Violette [1] is a single stage 
model that does not consider the oxygen pressure 
influence. Agarwal et al [3] suggested a homogeneous 
kinetic model considering the presence of the oxygen in 
the reaction media. 

The goal of this study was to obtain models 
adapted from previous published kinetics equations by 
using classical optimization methods (like: BFGS, DFP, 
Steepest Descent, Gauss Newton, Levenberg-Marquardt 
and Simplex) to represent the industrial process. After 
identifying the kinetic parameters as the Arrhenius 


frequency factor, activation energy and the exponents of 
the chemicals concentrations involved in the reaction, the 
models were applied in the commercial simulator 
CADSIM Plus.The software represents the process with 
mass and energy balances being capable to predict the 
kappa number after the oxygen delignification reactor. 

II. KINETIC MODELS 

Several authors studied and proposed kinetic 
models to the oxygen delignification reaction rate 
considering the mass transfer phenomena. Usually the 
reaction is described into two phases, the initial is fast and 
the final is slow or nonexistent lignin reaction. The fast 
reaction represents the alkaline extraction of the soluble 
lignin instead of the oxygen reactions [1]. This can be 
represented by a mathematical model composed of two 
parallel first order ordinary differential equations related to 
the lignin [4] 

The system can also be described by one equation 
representing the kappa number degradation in a potential 
form, usually having high reaction order related to the 
kappa number By using a simple reaction rate, the slow 
lignin degradation course during the final stage of the 
delignification can be mathematically described by the 
high reaction order. More slow the reaction, higher the 
exponent will be [4]. 

Agarwal et al [3] suggest representing the kinetic 
behavior as a series of parallel first order reactions 
occurring simultaneously and also consider the possibility 
of the nonexistent lignin presence characterizing the final 
slow reaction. The model is given by a one stage equation 
with a high reaction order related to the kappa number: 
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— = k-K 1J [OH ] 
dt 


r - - 1 0,92 


0,53 


( 1 ) 


1 07 2 

k = 2,368- 10 6 exp(-—) (2) 

RT 

The authors believe that there is some kind of 
delay in the process due to the series of parallel first order 
reactions that leads to a high order equation. The 
parameters identification based on industrial data is easy to 
be done using this model because it is represented by the 
global kappa number, not being necessary to characterize 
the fast and slow reaction lignin. This model is also 
efficient to apply into process simulators to evaluate the 
kappa number variations caused by the oxygen because the 
kinetic model considers the presence of this chemical in 
the reaction media. 

The objective of Violette’s [1] study was to 
improve the selectivity of the delignification reaction using 
polymeric additives. The author evaluated the effects of 
these additives in the selectivity and reaction kinetic rate. 
Violette [1] suggested that oxygen radicals may be trapped 
with relatively low amounts of additive radical scavengers 
if they are concentrated at the cellulosic surfaces. To 
confirm his hypotheses, the researcher determined the 
change in the kappa number given by the alkaline 
extraction by measuring the kappa variations in a reaction 
media without oxygen. 

The author tested seven models to evaluate the 
alkali consumption and he observed that the decrease in 
the alkali consumption was proportional to the kappa 
number decrease and independent of the process 
operational conditions: 

A[NaOH] = 0,168g / L • AK + 0,2g / L (3) 

The model to predict the kappa number is 
represented by equation 4. As Violette [1] used nitrogen 
instead of oxygen in his work, the model is represented by 
a corrected kappa number. As this study is based on 
industrial data the correction is not applicable. 

- — = k-K xn -[NaOH] 05SS (4) 

dt 

i a a ,^3 -7140. 

k = 4,4 • 10 3 • exp(-) (5) 

RT 

The parameter identification of Violette’s model 
is also simple to apply regarding the use of the global 
kappa number and a linear model of alkali consumption. 
However, as the oxygen is not considered in the kinetic 
equation rate, this model cannot be represented in the 


simulator as this does not show the kappa-oxygen 
interaction. 

III. PULP AND PAPER SIMULATORS 

The process simulation is a useful tool to reduce 
operational costs and improve the efficiency of existent 
mills, it is also essential to project new systems and make 
the start-ups faster and more rentable. Testing different 
solutions in the simulator helps to indentify the potential 
problems and change control strategy. 

During the 80’s, the computational system started 
to be trustable enough to simulate processes. The Fortran 
language gave the opportunity to realize complex 
calculations quickly. The solution was given in the 
stationary state, by mass and pressure flows, valve opening 
and tank levels. Then the process parameters were 
calculated, as concentration, reaction rate and temperature. 
In this same period, Honeywell acquired the software 
SCADA (Supervisory Control and Data Acquisition) 
which included solutions to differential algebraic equations 
and interaction with the control system. This resulted in 
engineer tools as the MASSBAL [5]. 

Since the early 90’s, the necessity of dynamic 
simulations appeared. Nowadays the most used pulp and 
paper simulators are the WinGEMS from Mesto 
Automation, IDEAS from Andritz and CADSIM Plus from 
Aurel Systems. 

CADSIM Plus allows a CAD draw to give the 
data source to the process simulation. It supports any 
drawing complexity, being especially adequate to P&ID 
diagrams. The dynamic simulation can be used to optimize 
process, to plan operations, to test DCS (Distributed 
Control System) control system and to develop control 
strategies. The software uses de ‘C’ code providing low 
execution time in commons computers. The user can also 
create particular modules of simulation called DLM 
(Dynamically Linked Modules). The streams of the 
process can indicate the fiber fraction and specific items as 
metallic ions, viscosity, brightness and kappa number [6]. 

The CADSIM Plus can use the function COM 
(Common Object Model) or DDE (Dynamic Data Change) 
to communicate with other applicative from Windows, as 
Microsoft Excel or Visual Basic. Once the communication 
is done, the software can send or receive data. The 
software can also Interact with distributes control systems 
(DCS) by the protocol OPC (OLE to process control) [6]. 
For example, in the Klabin Monte Alegre’s case, the 
communication with the PI software is possible. 

This simulator uses modules to represent the 
process. Each module is a mathematical representation of 
the process streams and equipment units. The software has 
a standard library that includes process equipments, 
control, integration, mathematical and logical modules. 
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There are optional libraries of fiber, utilities, 
hydrocarbons, data reconciliation and connection to other 
software and process controllers. The fiber library has 
modules for chip refinery, ion change, pressurized filters, 
pulp and reject refineries, washers and digesters washing 
zone. 


IV. DATA ACQUISITION AND MODELING 
DEVELOPMENT 

The data was collected from Klabin Papeis Monte 
Alegre industry The period was sufficient to recreate the 
process dynamics and generate the necessary amount of 
data to do trustable analyses. Statistical treatments were 
done to eliminate errors caused be production stops and 
variations, online analyzers wrong registered values and 
instruments out of operation. 

A total of 470 process variables patterns of pulp 
volumetric flow entering the stage, pressure and mass flow 
of the oxygen, reactor temperature, pulp consistency and 
kappa number from the Digester and inlet of the Bleach 
Plant were used to do the modeling and the process 
simulation. 

The kinetic parameters were adjusted based on the 
models given by Agarwal et al [3] and Violette [2]. The 
model to represent the alkali consumption is a linear model 
proportional do the kappa number decrease similar to the 
one proposed by Violette [2] given by equation 6. 

d[OH] dK 

-= a - (o) 

dt dt 

The first order differential equation system 
composed by the kinetic equation rate (equation 7 or 8 
referring to Violette and Agarwal’s et al model 
respectively) and the alkali consumption model was solved 
with direct integration with the fourth order Runge-Kutta 
method with variable step on the average residence time of 
the reactor. The parameter fitting was done using the 
classical optimization routines: BFGS, DFP, Steepest 
Descent, Gauss Newton, Levenberg-Marquardt and 
Simplex 


- —— = k • K a -[NaOH] b 
dt 

- — = k-K a [OH~] b P a c 
dt ° 2 


(7) 

( 8 ) 


To find the best equation to simulate the process 
the objective function was evaluated by the normalized 
quadratic error given by equation 9. A minimum value of 
the quadratic error implies that there is no need of a new 
adjusts of the parameters. 


W = -^-= (9) 

N-M 

The criteria adopted in the optimization routines 
of Matlab were a normalized error related to the calculated 
variable of 10-3 or a maximum number of iterations of 
750, what happened first. In all cases the test conjunct was 
33% form the original data to validate the models. 

The process simulation was created be drawing 
the industrial flowsheet with two main streams. One 
representing water, chemicals, organic and inorganic 
dissolved solids, fiber and lignin, the other is related to the 
steam line used in the process. 

The pulp washing was characterized by the 
displacement ratio and dilution factors particularized to 
each type of washer. The reactor simulation mode was 
done by mass balance based on the kinetic equation rates 
obtained after the optimization step. The equation rate is 
an input to the software that considers perfect mixture in 
the reaction vase. The type of the reactor chosen to 
represent the industrial process was an upflow tubular 
reactor. The tanks were modeled to represent a perfect 
mixture behavior. 

The software uses typical PID controllers to 
simulate de process control. The PID control varies the 
controlled variable until the measured variable reaches the 
set point which is a input to the simulation. CADSIM Plus 
tunes automatically the controllers but if the simulation 
shows oscillatory behavior or takes too long time to reach 
the set point it is possible to change the parameters values 
in order to avoid these disturbances. 

First, equation 7 with the adjusted parameters was 
used as the kinetic equation rate, followed by equation 8 
simulation. During the simulation, DDE communication 
protocols were created to make interactions between 
CADSIM Plus and Microsoft Excel. The simulator acts as 
data client or server, working as client the simulator 
imports the data to be utilized as inputs in the simulation. 
When working as server, the simulator exports the 
simulated data to the Microsoft Excel spreadsheets. 
Having the real process values and the simulated data, the 
relative error of the models was obtained. 

V. RESULTS 

To facilitate the comprehension of this work the 
models are nominated as Model 1, that is similar to 
Violette’s model and represents de first delignification 
line. For this same process line, Model 2 represent the one 
based on Agawrall et al study. Regarding the second 
delignification line, Model 3 represents the one modeled 
based on Violette’s proposition and Model 4 on 
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equations 10 and 11. The adjusted parameters were given 
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obtained by Rubini but the exponents of the chemicals 
concentrations changed. The difference between this and 
Rubini’s work is the raw material utilized in the 
delignification process. 
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- — = k- K 4 ’ 620 • [NaOH]~ 0 029 (10) 

dt 

-725-10 3 

fc = 11610-exp(---) (11) 

RT 

Figure 2 (A) and (B) illustrate the simulated 
results and the real industrial values. The first figure shows 
the values across the data quantity. The second one is the 
plotting of the calculated data by the real values, more the 
curve is close to the 45 line better the modeling prediction 
is. 



(A) 



(B) 

Fig. 2: Kappa number comparison between simulated and real values for Model 1. 


The simulated values follow the real values 
tendency even not considering the oxygen presence in the 
reaction media. Only for high kappa numbers the model 
was not accurate, however these are not usual operational 


values as they seldom occurs in the process. It can be 
concluded also that there is oxygen in excess because 
Violette’s model describes well the data pattern. 
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The absolute error (difference between real and 
calculated values) is + 2,5 kappa number and the relative 
error is 20%, the most part is due to process variations. 
The model is able to represent accurately normalized 
kappa numbers between 1,3 and 1,7. 

The kappa number decrease and the alkali 
consumption trough the reactor are illustrated on Figure 4. 
All optimization routines resulted in values close to 6 to 
the proportional constant that describes these profile, as 
given by equation 12, same result obtained by Rubini [2]. 

d[OH ] ^ dK 

-= 6 - ( 12 ) 

dt dt 

With this information, it is possible to evaluate 
the alkali amount consumed in the reaction and the 
possibility to reduce its concentration or flow in the 
process. These results are also normalized values. 




Fig.3: Alkali consumption and kappa number decrease 
profiles for Model 1. 

These profiles allow establishing quality limits to 
the reaction, i.e., find the optimum sodium hydroxide 
concentration to get the desired kappa number within the 
specifications [2]. The results show that the alkali 
consumption profile has similar behavior to the one 
propose by Violette [1]. 

Regarding Model 2, the optimization routines 
that gave the best results were Gauss Newton and 
Levenberg Marquardt. By choice, the results below are 
given by the Gauss Newton method. Equation 13 is the 
kinetic equation for this model, the kinetic constant is also 
represented by equation 11 because the results are the 
same as Model 1. The alkali consumption is the same as 
Model 1. 

- — = k- K 3330 -[NaOH] 0 ’ 151 (13) 

dt 


The absolute error is into the interval of + 2,5 
kappa number and the relative error is 15%, proving the 
model accuracy. Figure 4 (A) and (B) illustrate the 
simulated results and the real industrial values. 
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(A) 



(B) 

Fig. 4: Kappa number comparison between simulated and real values for Model 2. 


= 6,706.^ 

dt dt 


(14) 


The simulated values follow the real process 
values mainly in the range between 1,2 and 1,6 normalized 
kappa numbers, even not considering the oxygen presence 
in the reaction media. Again, it can be seen that the model 
is capable of predicting the process values as shown by 
Figure 4 (A). 

Model 3 was obtained by the Steepest Descent 
optimization method, having an relative error of 20% and 
an absolute error of + 2,0 kappa number. Equation 14 
represents the alkali consumption and equations 15 and 16 
show the kinetic rate. 


— = k • K 9 ’ 231 • [OH-] 5 * 39 • P Q 41 ’ 986 (15) 

dt ° 2 

, „„ , -107,20. 

k = 2,3 -10 6 -exp(-) (16) 

RT 
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Figure 5 (A) and (B) illustrate this model is able 1,2 and 1,6, only high kappa numbers were not well 

to predict normalized kappa number values in the range of represented but they are not usual in the process. 



(A) 



(B) 

Fig. 5: Kappa number comparison between simulated and real values for Model 3. 


Model 4 was better represented by the BFGS 
method as this was the one with the lowest error conjunct 
data. Equations 17, 18 and 19 represent the kinetic rate and 
the alkali consumption model respectively with accuracy 
of ± 2,0 kappa number and 15% of relative error. 


— = k • K 1MA • [OH~] 2,694 • P Q -3 ’ 343 (17) 

dt ° 2 

-72 5-10 3 

k = 11610exp( ---) (18) 

RT 
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^^ = 6 , 612 - — 
dt dt 


As shown by the figures below, this model has a 
(1^) good performance between the range of 1,2 and 1,6 

normalized kappa number. The alkali consumption profile 
is similar to the one from Model 1. 



(A) 



(B) 

Fig. 6: Kappa number comparison between simulated and real values for Model 4. 
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Only Models 3 and 4 gave coherent result when 
applied to the CADSIM Plus simulator, also these models 
had short simulation time. 

In the first simulation, the software acted as client 
receiving the inputs, as chemical flows and concentrations 



from excel. It presented an average absolute error of 1,3 
kappa number after the delignification reactor. Figure 7 is 
a tendency chart of the real and simulated kappa number. 


Fig. 7: Caparison between real and simulated kappa number values simulated on CADSIM Plus with Model 3. 


The figure above shows that this model is able to 
follow the process tendency, excepting for some points, 
probably from transient periods in the operation or 
conditions that were not used in the data fitting. Even with 
the data treatment to assure stationary simulation, some 
data should be characteristic of process variations. This is 
because this is a new operation in the mill, resulting in 
some instability 

With CADSIM Plus it was possible to evaluate de 
chemical consumption inside the delignification reactor, 


simulating the mass balance to the oxygen with mass flow 
into the process and the residual oxygen after the reaction, 
the same was done referring to alkali. The result showed 
oxygen in excess while the alkali was all consumed in the 
reaction corroborating with the results obtained from the 
alkali consumption modeling. 

Model 4 simulation gave an average absolute 
error of 1,2 kappa number. Figure 8 shows the tendency 
chart of the real and simulated kappa number. 



Fig. 8: Caparison between real and simulated kappa number values simulated on CADSIM Plus with Model 4. 
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Again, this model is able to follow the behavior of 
the process excepting for unusual data. When analyzing 
the mass balance for oxygen and alkali, it ca be seen that 
this line also operates with excess oxygen while practically 
all the alkali is consumed. 

In the second type of simulation for both models, 
the software worked as server, with fixed usual process 
conditions as inputs, providing the analyses of the kappa 
number behavior when changing operational conditions. 

VI. CONCLUSIONS 

To represent the kinetic of the chemical reactions 
two adapted phenomenological models based on literature 
study were used. Both models were able to predict the 
kappa number after the delignification reactor and it was 
possible to obtain profiles for the kappa number decrease 
and the alkali consumption inside the reactor. However 
only the model that considered the oxygen presence in the 
reaction media was able to be implemented in the 
CADSIM Plus simulator generating reliable results. 

The vantage of the phenomenological modeling is 
that, once the models are validated, they become a 
methodology potentially useful to investigate the behavior 
of the reaction in many operational conditions. 

The models showed an average relative error of 
13% to one delignification line and 10% to the other. The 
models are applicable to predict normalized kappa 
numbers in the range of 1,2 to 1,6 kappa number units. 
Rubini’s model presented an average relative error of 
8,5%, these values are in agreement with the errors given 
by the literature, about 8 to 20% [2]. 

This shows that the models obtained in this study 
have good accuracy to represent the real process as they 
are in the same level of the errors obtained by other 
authors in the literature. The deviations were slightly 
higher from another works because this is a new process in 
the mill operation what results in process variations and 
non stationary operation. Even with cautiously analyses to 
assure the stationary state, some data still can have errors 
from instrument calibration, data acquisition, or signal 
processing for PI software. 
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